/**
 * Mathematic : Definition de la loi Gamma
 *
 * @author maria ladoue
 * @version 0.1, 30/08/02 
 */
package br.unb.statistic;

import VisualNumerics.math.Statistics;
import br.unb.utils.Resource;

public class Mathematic implements ComputerArithmetic {

    //
    // Relative to polynom evaluation
    //
    /**  
     *    Evaluate polynomial
     *
     * @param x : value 
     * @param coeff : coefficient of corresponding degree
     * @param N : degree
     *
     * y = polEval( x, coef, N );
     *
     *
     * DESCRIPTION:
     *
     * Evaluates polynomial of degree N:
     *
     * y(x)  =  C  + C x + C x^2  +...+ C x^N
     *           0    1     2            N
     *
     * Coefficients are stored in reverse order:
     *
     * coef[0] = C  , ..., coef[N] = C  .
     *            N                   0
     *
     *  The function p1evl() assumes that coef[N] = 1.0 and is
     * omitted from the array.  Its calling arguments are
     * otherwise the same as polevl().
     *
     * SPEED:
     *
     * In the interest of speed, there are no checks for out
     * of bounds arithmetic.  This routine is used by most of
     * the functions in the library.  Depending on available
     * equipment features, the user may wish to rewrite the
     * program in microcode or assembly language.
     *
     * Cephes Math Library Release 2.1:  December, 1988
     * Copyright 1984, 1987, 1988 by Stephen L. Moshier
     * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
     */
    public static double polEval(double x, double[] coef, int N) {

        double ans = 0;

        if (N >= 0) {
            for (int i = N; i >= 0; i--) {
                // keep multiplying by x and adding : a0 + x( a1 + x (a2 ...  )))
                ans = ans * x + coef[i];
            }
            return ans;  // to perform
        } else {
            throw new IllegalArgumentException(Resource.getString("MATHEMATIC.P1EVL_MUST_BE_>_1"));
        }
    } //end of polEval

    /**                                        
     * Evaluate polynomial when coefficient of x^N is 1.0 => Cn=1.
     * Otherwise same as polEval.
     */
    public static double pol1Eval(double x, double[] coef, int N) throws IllegalArgumentException {

        double ans = 1;

        if (N >= 1) {
            for (int i = N - 1; i >= 0; i--) {
                // keep multiplying by x and adding : a0 + x( a1 + x (a2 ...  )))
                ans = ans * x + coef[i];
            }
            return ans;
        } else {
            throw new IllegalArgumentException(Resource.getString("MATHEMATIC.P1EVL_MUST_BE_>_1"));
        }
    } //end of pol1Eval

    //
    //  Relative to : Normal Function
    //
    /**
     *      Inverse of Normal distribution function
     *
     * @param x double
     * @param y double
     *
     * x = invNormal( y );
     *
     * DESCRIPTION:
     *
     * Returns the argument, x, for which the area under the
     * Gaussian probability density function (integrated from
     * minus infinity to x) is equal to y.
     *
     *
     * For small arguments 0 < y < exp(-2), the program computes
     * z = sqrt( -2.0 * log(y) );  then the approximation is
     * x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
     * There are two rational functions P/Q, one for 0 < y < exp(-32)
     * and the other for y up to exp(-2).  For larger arguments,
     * w = y - 0.5, and  x/sqrt(2pi) = w + w^3 R(w^2)/S(w^2)).
     *
     *
     * ACCURACY:
     *
     *                      Relative error:
     * arithmetic   domain        # trials      peak         rms
     *    DEC      0.125, 1         5500       9.5e-17     2.1e-17
     *    DEC      6e-39, 0.135     3500       5.7e-17     1.3e-17
     *    IEEE     0.125, 1        20000       7.2e-16     1.3e-16
     *    IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
     *
     *
     * ERROR MESSAGES:
     *
     *   message         condition    value returned
     * invNormal domain       x <= 0        -MAXNUM
     * invNormal domain       x >= 1         MAXNUM
     *
     * Cephes Math Library Release 2.8:  June, 2000
     * Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
     */
    public static double invNormal(double y0) {

        return Statistics.inverseNormalCdf(y0);

        /**
        double x, y, z, y2, x0, x1;
        int code;

        if( y0 <= 0.0 ){
        //System.out.println("Mathematic.invNormal : DOMAIN ERROR");
        return  -MAXNUM;
        }
        if( y0 >= 1.0 ) {
        //System.out.println("Mathematic.invNormal : DOMAIN ERROR");
        return MAXNUM ;
        }
        code = 1;
        y = y0;
        if( y > 1.0 - 0.13533528323661269189 ) { // 0.135... = exp(-2)
        y = 1.0 - y;
        code = 0;
        }

        if( y > 0.13533528323661269189 ) {
        y = y - 0.5;
        y2 = y * y;
        x = y + y * (y2 * polEval( y2, P0, 4) / pol1Eval( y2, Q0, 8 ));
        x = x * s2pi;
        return x;
        }

        x = Math.sqrt( -2.0 * Math.log(y) );
        x0 = x - Math.log(x)/x;

        z = 1.0/x;
        if( x < 8.0 ) {  // y > exp(-32) = 1.2664165549e-14
        x1 = z * polEval( z, P1, 8 ) / pol1Eval( z, Q1, 8 );
        }
        else {
        x1 = z * polEval( z, P2, 8 ) / pol1Eval( z, Q2, 8 );
        }
        x = x0 - x1;
        if( code != 0 ) { x = -x; }

        return x ;

         */
    } // end of invNormal

    //
    //  Relative to : Gamma Function
    //
    /**                                                     
     *      Inverse of complemented imcomplete gamma integral
     *
     * @param a
     * @param y0
     *
     * DESCRIPTION:
     *
     * Given p, the function finds x such that
     *
     *  incompleteGammaComplemented( a, x ) = p.
     *
     * Starting with the approximate value
     * 
     *  x = a t^3
     *
     *  where :
     *
     *    t = 1 - d - invNormal(p) sqrt(d)   and  d = 1/9a,
     *
     * It performs up to 10 Newton iterations to find the
     * root of incompleteGammaComplemented(a,x) - p = 0.
     *
     * ACCURACY:
     *
     * Tested at random a, p in the intervals indicated.
     *
     *                a        p                      Relative error:
     * arithmetic   domain   domain     # trials      peak         rms
     *    IEEE     0.5,100   0,0.5       100000       1.0e-14     1.7e-15
     *    IEEE     0.01,0.5  0,0.5       100000       9.0e-14     3.4e-15
     *    IEEE    0.5,10000  0,0.5        20000       2.3e-13     3.8e-14
     *
     * Cephes Math Library Release 2.8:  June, 2000
     * Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
     */
    public static double incompleteGammaComplementedInverse(double a, double y0) {

        double x0, x1, x, yl, yh, y, d, lgm, dithresh;
        int dir;
        boolean intervalHalving = false, done = false;

        // bound the solution
        x0 = MAXNUM;
        yl = 0;
        x1 = 0;
        yh = 1.0;
        dithresh = 5.0 * MACHEP;

        // approximation to inverse function
        d = 1.0 / (9.0 * a);
        y = (1.0 - d - invNormal(y0) * Math.sqrt(d));
        x = a * y * y * y;

        lgm = lnGamma(a);

        for (int i = 0; i < 10; i++) {
            if (x > x0 || x < x1) {
                intervalHalving = true;
                break;
            }

            y = incompleteGammaComplemented(a, x);
            if (y < yl || y > yh) {
                intervalHalving = true;
                break;
            }
            if (y < y0) {
                x0 = x;
                yl = y;
            } else {
                x1 = x;
                yh = y;
            }
            // compute the derivative of the function at this point
            d = (a - 1.0) * Math.log(x) - x - lgm;
            if (d < -MAXLOG) {
                intervalHalving = true;
                break;
            }
            d = -Math.exp(d);
            // compute the step to the next approximation of x
            d = (y - y0) / d;
            if (Math.abs(d / x) < MACHEP) {
                done = true;
                break;
            }

            x = x - d;
            return x;  //debug maria 12/09/02
        }

        //Resort to interval halving if Newton iteration did not converge.

        if (intervalHalving) {

            d = 0.0625;
            if (x0 == MAXNUM) {

                if (x <= 0.0) {
                    x = 1.0;
                }

                while (x0 == MAXNUM) {
                    x = (1.0 + d) * x;
                    y = incompleteGammaComplemented(a, x);
                    if (y < y0) {
                        x0 = x;
                        yl = y;
                        break;
                    }
                    d = d + d;
                }
            }

            d = 0.5;
            dir = 0;

            for (int i = 0; i < 400; i++) {
                x = x1 + d * (x0 - x1);
                y = incompleteGammaComplemented(a, x);

                lgm = (x0 - x1) / (x1 + x0);
                if (Math.abs(lgm) < dithresh) {
                    break;
                }

                lgm = (y - y0) / y0;
                if (Math.abs(lgm) < dithresh) {
                    break;
                }

                if (x <= 0.0) {
                    break;
                }

                if (y >= y0) {
                    x1 = x;
                    yh = y;
                    if (dir < 0) {
                        dir = 0;
                        d = 0.5;
                    } else if (dir > 1) {
                        d = 0.5 * d + 0.5;
                    } else {
                        d = (y0 - yl) / (yh - yl);
                    }
                    dir += 1;
                } else {
                    x0 = x;
                    yl = y;
                    if (dir > 0) {
                        dir = 0;
                        d = 0.5;
                    } else if (dir < -1) {
                        d = 0.5 * d;
                    } else {
                        d = (y0 - yl) / (yh - yl);
                    }
                    dir -= 1;
                }
            }  // end of for

            if (x == 0.0) {
                throw new ArithmeticException(Resource.getString("MATHEMATIC_INCOMPLETE_GAMMA_COMPLEMENTED_INVERSE_OVERFLOW"));
            }

        }  // end of intervalHalving

        if (done) {
            return x;
        }
        return x;//debug maria 12/09/02

    }// end of incompleteGammaComplementedInverse
    /**
     *
     *      Incomplete gamma integral
     *
     *
     *
     * SYNOPSIS:
     *
     * double a, x, y, igam();
     *
     * y = igam( a, x );
     *
     * DESCRIPTION:
     *
     * The function is defined by
     *
     *                           x
     *                            -
     *                   1       | |  -t  a-1
     *  igam(a,x)  =   -----     |   e   t   dt.
     *                  -      | |
     *                 | (a)    -
     *                           0
     *
     *
     * In this implementation both arguments must be positive.
     * The integral is evaluated by either a power series or
     * continued fraction expansion, depending on the relative
     * values of a and x.
     *
     * ACCURACY:
     *
     *                      Relative error:
     * arithmetic   domain     # trials      peak         rms
     *    IEEE      0,30       200000       3.6e-14     2.9e-15
     *    IEEE      0,100      300000       9.9e-14     1.5e-14
     */
    /**
     *
     *      Complemented incomplete gamma integral
     *
     *
     *
     * SYNOPSIS:
     *
     * double a, x, y, igamc();
     *
     * y = igamc( a, x );
     *
     * DESCRIPTION:
     *
     * The function is defined by
     *
     *
     *  igamc(a,x)   =   1 - igam(a,x)
     *
     *                            inf.
     *                              -
     *                     1       | |  -t  a-1
     *               =   -----     |   e   t   dt.
     *                    -      | |
     *                   | (a)    -
     *                             x
     *
     *
     * In this implementation both arguments must be positive.
     * The integral is evaluated by either a power series or
     * continued fraction expansion, depending on the relative
     * values of a and x.
     *
     * ACCURACY:
     *
     * Tested at random a, x.
     *                a         x                      Relative error:
     * arithmetic   domain   domain     # trials      peak         rms
     *    IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
     *    IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15
     *
     * Cephes Math Library Release 2.8:  June, 2000
     * Copyright 1985, 1987, 2000 by Stephen L. Moshier
     */
    final static double big = 4.503599627370496e15;
    final static double biginv = 2.22044604925031308085e-16;

    public static double incompleteGammaComplemented(double a, double x) {
        double ans, ax, c, yc, r, t, y, z;
        double pk, pkm1, pkm2, qk, qkm1, qkm2;

        if ((x <= 0) || (a <= 0)) {
            return (1.0);
        }

        if ((x < 1.0) || (x < a)) {
            return (1.0 - incompleteGamma(a, x));
        }

        ax = a * Math.log(x) - x - lnGamma(a);
        if (ax < -MAXLOG) {
            //mtherr( "igamc", UNDERFLOW );
            return (0.0);
        }
        ax = Math.exp(ax);

        // continued fraction
        y = 1.0 - a;
        z = x + y + 1.0;
        c = 0.0;
        pkm2 = 1.0;
        qkm2 = x;
        pkm1 = x + 1.0;
        qkm1 = z * x;
        ans = pkm1 / qkm1;

        do {
            c += 1.0;
            y += 1.0;
            z += 2.0;
            yc = y * c;
            pk = pkm1 * z - pkm2 * yc;
            qk = qkm1 * z - qkm2 * yc;
            if (qk != 0) {
                r = pk / qk;
                t = Math.abs((ans - r) / r);
                ans = r;
            } else {
                t = 1.0;
            }
            pkm2 = pkm1;
            pkm1 = pk;
            qkm2 = qkm1;
            qkm1 = qk;
            if (Math.abs(pk) > big) {
                pkm2 *= biginv;
                pkm1 *= biginv;
                qkm2 *= biginv;
                qkm1 *= biginv;
            }
        } while (t > MACHEP);

        return (ans * ax);
    }

    /**
     *left tail of incomplete gamma function:
     *
     *          inf.      k
     *   a  -x   -       x
     *  x  e     >   ----------
     *           -     -
     *          k=0   | (a+k+1)
     *
     */
    public static double incompleteGamma(double a, double x) {
        double ans, ax, c, r;

        if ((x <= 0) || (a <= 0)) {
            return (0.0);
        }

        if ((x > 1.0) && (x > a)) {
            return (1.0 - incompleteGammaComplemented(a, x));
        }

        // Compute  x^a * exp(-x) / gamma(a)
        ax = a * Math.log(x) - x - lnGamma(a);
        if (ax < -MAXLOG) {
            // mtherr( "igam", UNDERFLOW );
            return (0.0);
        }
        ax = Math.exp(ax);

        /* power series */
        r = a;
        c = 1.0;
        ans = 1.0;

        do {
            r += 1.0;
            c *= x / r;
            ans += c;
        } while (c / ans > MACHEP);

        return (ans * ax / a);
    }

    //
    //  Methode permettant de calculer la fonction Gamma Incomplete
    //     P(a,x)=1/Gamma(a)*Int[exp(-t)*t^(a-1)dt,0,x] pour a>0
    //    Q(a,x)~1-P(a,x)~1/Gamma(a)*Int[exp(-t)*t^(a-1)dt,x,infinity] pour a>0
    //  Methode recuperer de NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING
    //
    /**
     * Incomplete Gamma Integral
     */
    /**   public static double incompleteGamma (double a,double x) throws IllegalArgumentException {
    double gln;
    if(x < 0d || a<=0d) {
    throw new IllegalArgumentException("GammaDistribution.iGammaP : Invalid arguments");
    }
    else if(x < a+1d) {    // Use the series representation
    return gser(a,x);
    }
    else {
    return 1d - gcf(a,x);
    }
    }
     */
    /**
     * Complemented Incomplete Gamma Integral
     */
    /**
    public static double incompleteGammaComplemented (double a,double x) {
    double gamser,gammcf,gln;
    if(x < 0d || a <= 0d) {  // Use the series representation
    throw new IllegalArgumentException("GammaDistribution.iGammaP : Invalid arguments");
    }
    else if (x < (a+1d)) {    // Use the series representation
    return 1-gser(a,x);
    }
    else {
    return gcf(a,x);
    }

    }// end of 
     */
    /**
     *
     * Returns the logarithm of the Gamma function of a double.
     *
     * @param xx  A double value
     *
     *  Returns the value of ln[Gamma(xx)] for xx>0 
     *  Methodes takes from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING
     *  If x is a negative integer, the result is NaN.
     */
    public static double lnGamma(double xx) throws IllegalArgumentException {
        if (xx > 0d) {
            double x, y, tmp, ser;
            final double[] cof = {76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179 * Math.exp(-2), -0.5395239384953 * Math.exp(-5)};

            y = x = xx;
            tmp = x + 5.5;
            tmp -= (x + 0.5) * Math.log(tmp);
            ser = 1.000000000190015;
            for (int j = 0; j <= 5; j++) {
                ser += cof[j] / ++y;
            }
            return -tmp + Math.log(2.5066282746310005 * ser / x);
        } else {
            return Double.NaN;
            //throw new IllegalArgumentException("GammaDistribution.lnGamma : argument must be > 0.");
        }
    } // end of lnGamma
    // Maximum allowed numbre of iterations
    final static int ITMAX = 100;
    // Relative accuracy
    final static double EPS = 3.0 * Math.exp(-7);
    // Numbre near the smallest representable floating-point number
    final static double FPMIN = Double.MIN_VALUE;  // MIN_VALUE = 4.9E-324d dans l'algo 1.0*Math.exp(-30);//

    /**
     *  Returns the incomplete gamma function P(a,x) evaluated by its series representation as gamser 
     */
    public static double gser(double a, double x) throws IllegalArgumentException {

        double sum, del, ap;
        double gln = lnGamma(a);
        double gamser;

        if (x <= 0d) {
            if (x < 0d) {
                throw new IllegalArgumentException(Resource.getString("GAMMADISTRIBUTION.GSER_X_LESS_THAN_0"));
            }
            gamser = 0d;
            return gamser;
        } else {
            ap = a;
            del = sum = 1d / a;
            for (int n = 1; n <= ITMAX; n++) {
                ++ap;
                del *= x / ap;
                sum += del;
                if (StrictMath.abs(del) < StrictMath.abs(sum) * EPS) {
                    gamser = sum * Math.exp(-x + a * Math.log(x) - gln);
                    return gamser;
                }
            }
            throw new IllegalArgumentException(Resource.getString("GAMMADISTRIBUTION.GSER_A_TOO_LARGE,_ITMAX_TOO_SMALL."));
        }

    }// end of gser

    /**
     * Returns the incomplete gamma function Q(a,x) evaluated by its continued fraction representation as gammcf.
     */
    public static double gcf(double a, double x) throws IllegalArgumentException {

        double gammcf;
        double an, b, c, d, del, h;
        int i;
        double gln = lnGamma(a);

        // set up for evaluating continued fraction by modified Lentz's method.
        b = x + 1d - a;
        c = 1d / FPMIN;
        d = 1.0 / b;
        h = d;
        for (i = 1; i <= ITMAX; i++) {   // Iterate to convergence
            an = -i * (i - a);
            b += 2d;
            d = an * d + b;
            if (StrictMath.abs(d) < FPMIN) {
                d = FPMIN;
            }
            c = b + an / c;
            if (StrictMath.abs(c) < FPMIN) {
                c = FPMIN;
            }
            d = 1d / d;
            del = d * c;
            h *= del;
            if (StrictMath.abs(del - 1d) < EPS) {
                break;
            }
        }

        if (i > ITMAX) {
            throw new IllegalArgumentException(Resource.getString("GAMMADISTRIBUTION.GCF_A_TOO_LARGE,_ITMAX_TOO_SMALL."));
        }
        gammcf = Math.exp(-x + a * Math.log(x) - gln) * h;  //Put factors in front
        return gammcf;

    }  // end of gcf

    /**
     * Gamma function computed by Stirling's formula.
     * The polynomial STIR is valid for 33 <= x <= 172.
     */
    public static double stirlingFormula(double x) {
        double y, w, v;

        w = 1d / x;
        w = 1d + w * polEval(w, STIR, 4);
        y = Math.exp(x);
        if (x > MAXSTIR) { // Avoid overflow in Math.pow()
            v = Math.pow(x, 0.5 * x - 0.25);
            y = v * (v / y);
        } else {
            y = Math.pow(x, x - 0.5) / y;
        }
        y = SQTPI * y * w;
        return y;
    }

    /**
     * Method factorial : Factorial function: compute factorial of n
     *
     * @param n  it'a a long
     * @return n!
     */
    public static double factorial(double n) throws IllegalArgumentException {
        if (n >= 0) {
            double ans = 1;
            for (double i = 2; i <= n; i++) {
                ans = ans * i;
            }
            if (ans == 0.0) {
                throw new IllegalArgumentException(
                        Resource.getString("MATHEMATIC.FACTORIAL: ")
                        + n
                        + Resource.getString("IS_TOO_BIG_FOR_JAVA_FACTORIAL_COMPUTATION!_N_MAX_IS_170!"));
            } else {
                return ans;
            }
        } else {
            throw new IllegalArgumentException(Resource.getString("MATHEMATIC.FACTORIAL_N_0,1,...N"));
        }
    }

    /**
     * Method pgcd : calcule le PGCD de 2 entiers passes en arguments (on utilisera le fait que le PGCD de 2 nbs divise le reste de la divisionm de ces 2 nombres en iterant)
     * @param a
     * @param b
     * @return
     */
    public static int pgcd(int a, int b) {
        if (a < b) {
            return pgcd(b, a);
        }
        while (b != 0) {
            int r = a % b;
            a = b;
            b = r;
        }
        return a;
    }

    /**
     * Calcule le PPCM de deux entiers passes en arguments
     * @param a an int
     * @param b an int
     * @return an int
     */
    public static int ppcm(int a, int b) {
        return (a * b) / pgcd(a, b);
    }

    /**
     * Method ppcmComplex : to compute the MCM(Minimum Commun Multiple) of a table of int
     * 		Property: MCM(A,B,C) =  MCM( MCM(A,B),C)....
     * @param nbs
     * @return an int
     */
    public static int ppcmComplex(int[] nbs) {
        int ppcm = nbs[0];

        for (int i = 1; i < nbs.length; i++) {
            ppcm = Mathematic.ppcm(ppcm, nbs[i]);
        }
        return ppcm;
    }

    /**
     * Methode: Main pour tester la classe
     */
/*
    public static void main(String[] args) {

        System.out.println( Mathematic.incompleteGammaComplemented( 2,0.1 ));  //OK
	System.out.println( Mathematic.lnGamma(0.2)); //OK
	System.out.println( Mathematic.lnGamma(0.5)); //OK
	System.out.println( Mathematic.lnGamma(10));
	System.out.println( Mathematic.lnGamma(15));
	System.out.println( Mathematic.lnGamma(45));

	System.out.println( Sfun.logGamma(0.2)); //OK
	System.out.println( Sfun.logGamma(0.5)); //OK
	System.out.println( Sfun.logGamma(10));
	System.out.println( Sfun.logGamma(15));
	System.out.println( Sfun.logGamma(45));

        System.out.println( Mathematic.invNormal(0.2));
        System.out.println( Mathematic.invNormal(0.5));
        System.out.println( Mathematic.invNormal(0.7));
	System.out.println( Mathematic.factorial(0));
	System.out.println( Mathematic.factorial(100));

        System.out.println("PPCM: " + Mathematic.ppcm(2, 3));
        System.out.println("PPCM: " + Mathematic.ppcm(3, 5));
        System.out.println("PPCM: " + Mathematic.ppcm(2, 5));
	System.out.println( "PPCM: " + Mathematic.ppcm(250, 400) );

	System.out.println( "PGCD: " + Mathematic.pgcd(480, 250) );
	System.out.println( "PGCD: " + Mathematic.pgcd(480, 300) );
	System.out.println( "PGCD: " + Mathematic.pgcd(250, 300) );
	System.out.println( "PGCD: " + Mathematic.pgcd(250, 400) );
	
        int[] t = {2, 3, 5};
        System.out.println("PPCMComplex: " + Mathematic.ppcmComplex(t));
    }*/
}
